function [estPoint] = mullersMethod(func,threeApproxPoints,tol,maxIter)
%MULLERSMETHOD Function that determines the root of a linear function
%using Mullers Method, particularly for finding roots of polynomials.
%   This algorithm was adapted from Chapter 2.6 of Numerical Analysis (8th
%   edition, Burden & Faires)
%   
%   *Input Parameters*
%   func: String containing name of function to be evaluated
%   threeApproxPoints: Array of three approximate values of the zero of the
%                      function.
%   tol: Tolerance limit. 
%   maxIter: Maximum number of iterations
%
%   *Default options* 
%   tol = 10^-5
%   maxIter = 100
%
%   ***If no fixed point is found after maxIteration, the midpoint for 
%   maxIter approximation is returned.***
%
%   ***********************************************************************
%   Author: Mathieu Boudreau, BSc, MSc, PhD Candidate (BME)
%   Institute: Montreal Neurological Institute, McGill University
%   Contact: mathieu.boudreau2 (at) mail.mcgill.ca
%   Date: July 17th 2014
%   ***********************************************************************

%% Set default input conditions, if required.
%

if nargin < 4
    maxIter=100;
end

if nargin < 3
    tol=10^-5;
end


%% Initialize initial conditions
%
p0=threeApproxPoints(1);
p1=threeApproxPoints(2);
p2=threeApproxPoints(3);

h1=p1-p0;
h2=p2-p1;

del1=(feval(func,p1)-feval(func,p0))/h1;
del2=(feval(func,p2)-feval(func,p1))/h2;

d=(del2-del1)/(h2+h1);

initIter=3;


%% Run Muellers Algorithm
%

for iter=initIter:maxIter
   
    b=del2+h2*d;
    D=sqrt(b^2-4*feval(func,p2)*d);
    
    if sign(b)>0
        E=b+D;
    else
        E=b-d;
    end
    
    h=-2*feval(func,p2)/E;
    p=p2+h;
    
    if abs(h)<tol
        estPoint=p;
    	disp(['Algorithm terminated. Tolerance condition was reached after ' num2str(iter) ' iterations.'])
        return
    end
    
    p0=p1;
    p1=p2;
    p2=p;

    h1=p1-p0;
    h2=p2-p1;

    del1=(feval(func,p1)-feval(func,p0))/h1;
    del2=(feval(func,p2)-feval(func,p1))/h2;

    d=(del2-del1)/(h2+h1);
end

if iter==maxIter
    estPoint=p;
    disp('Function terminated after the maximum number iterations. No solution found within tolerance.')
end

end
